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\  ABSTRACT 

V 

A  decomposition  procedure  is  proposed  in  this  paper  for  solving  a 
class  of  large  scale  optimal  design  problems  for  perfectly  plastic  structures 
under  several  alternative  loading  conditions.  The  conventional  finite 
element  method  is  used  to  cast  the  problem  into  a  finite  dimensional  con¬ 
strained  nonlinear  programming  problem.  Structures  of  practically  meaningful 
size  and  complexity  tend  to  give  rise  to  a  large  number  of  variables  and 
constraints  in  the  corresponding  mathematical  model.  The  difficulty  is  that 
the  state-of-the-art  Mathematical  Programming  theory  does  not  provide 
reliable  and  efficient  ways  of  solving  large  scale  constrained  nonlinear 
programming  problems.  The  natural  idea  to  deal  with  the  large  scale 
structural  problem  is  to  somehow  decompose  the  problem  into  an  assembly  of 
small  size  problems  each  of  which  represents  an  analysis  of  the  behavior  of 
each  finite  element  under  a  single  loading  condition.  This  paper  proposes 
one  such  way  of  decomposition  based  on  the  duality  theory  and  a  recently 
developed  iterative  algorithm. 
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SIGNIFICANCE  AND  EXPLANATION 


Many  optimization  and  design  problems  arising  in  structural  analysis  are 
formulated  and  solved  as  nonlinear  programming  problems.  A  variety  of 
formulations  and  accompanying  solution  methods  have  been  proposed.  Structures 
with  practically  meaningful  size  and  complexity  tend  to  give  rise  to  a  large 
number  of  constraints  and  variables  in  the  corresponding  mathematical  model. 
This  tendency  is  particularly  prominent  when  one  uses  the  conventional  Finite 
Element  Method,  where  a  structure  is  discretized  into  a  large  number  of  small 
"finite"  elements.  A  major,  unresolved  difficulty  is  that  the  state-of-the-art 
nonlinear  programming  theory  does  not  provide  a  reliable,  efficient  way  of 
solving  a  programming  problem  with  a  large  number  of  nonlinear  constraints. 

This  paper  proposes  a  decomposition  procedure  for  a  limit  analysis  of  a  certain 
broad  class  of  (perfectly-plastic)  structures.  The  main  benefit  of  this 
procedure  is  that  it  allows  one  to  solve  (many)  small-size  nonlinear  program¬ 
ming  problems,  each  of  which  corresponds  to  an  individual  finite  element, 
instead  of  having  to  solve  the  entire  problem  involving  a  potentially  huge 
number  of  constraints  and  variables. 
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A  DECOMPOSITION  PROCEDURE  FOR  LARGE  SCALE 
OPTIMAL  PLASTIC  DESIGN  PROBLEMS 


Ikuyo  Kaneko*  and  Cu  Duong  Ha** 

1 .  Introduction 

In  this  paper  we  shall  consider  an  optimal  design  problem  of  perfectly 
plastic  structures  of  potentially  large  sizes  under  fairly  general  assumptions. 
Using  the  conventional  finite  element  method  the  structure  is  discretized  into 
an  assembly  of  elements,  and  the  problem  is  formulated  as  a  constrained  non¬ 
linear  programming  problem  of  finite  dimension.  A  number  of  ways  of  formulat¬ 
ing  the  optimal  design  problem  into  a  mathematical  programming  problem  have 
been  proposed.  The  unresolved  difficulty  lies  in  the  facts  that  (i)  the  pro¬ 
gramming  problems  arising  from  real-world  structural  problems  tend  to  have 
very  large  sizes  and  (ii)  the  state-of-the-art  Mathematical  Programming  theory 
does  not  provide  reliable  and  effective  algorithms  to  solve  large  size  con¬ 
strained  nonlinear  problems. 

Using  the  finite  element  method  it  is  inevitable  that  the  number  of 
elements  becomes  large  (hundreds  and  maybe  thousands)  if  we  want  to  handle 
engineering  structures  of  practically  meaningful  size  and  complexity  with 
reasonable  accuracy.  In  most  cases,  each  finite  element  is  associated  with  a 
relatively  small  number  of  variables  (generalized  stresses,  design  variables, 
etc.)  and  constraints  (yield  conditions,  etc.),  but  because  of  obvious  inter¬ 
actions  among  elements  (e.g.,  the  overall  equilibrium  condition)  we  need  to 
consider  all  the  finite  elements  simultaneously,  resulting  in  a  potentially 
huge  number  of  variables  and  constraints. 
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Much  research  is  being  done  in  the  field  of  Large  Scale  Mathemat¬ 
ical  Programming  but  at  this  point,  virtually  no  reliable,  well-tested 

algorithms  are  available  to  solve  optimization  problems  with  nonlinear 
constraints  having  more  than,  say,  a  hundred  variables  and  constraints. 

It  must  be  pointed  out  that  prospect  is  much  brighter  for  unconstrained 
minimization  of  a  differentiable  convex  function  (or  maximization  of 
differentiable  concave  function)  especially  when  the  explicit  form  of 
derivatives  are  known.  In  such  a  case,  some  of  the  most  reliable  and 
fast  algorithms  (such  as  quasi-Newton  methods)  can  be  used  to  handle 
large  size  problems  in  a  reasonably  effective  manner.  Of  course,  however, 
the  optimal  design  problems  of  structures  do  have  nonlinear  constraints 
and  thus  these  unconstrained  techniques  can  not  be  applied,  at  least 
directly. 

These  considerations  lead  to  the  following  intuitive  idea:  Instead 
of  dealing  with  one  large  nonlinear  program  corresponding  to  the  entire 
structure,  can  one  decompose  the  problem  into  small  size  programs  each 
of  which  corresponds  to  the  analysis  of  the  behavior  of  the  individual 
finite  elements  under  a  single  loading  condition  and  somehow  combine  the 
results  of  these  analyses  to  produce  the  solution  to  the  entire  structure? 
We  believe  that  this  is  a  natural  and  effective  approach  to  handle  large 
scale  structural  optimization  problems;  and  in  this  paper  we  shall  propose 
one  such  approach. 

The  following  is  an  informal  outline  of  the  proposed  decomposition 
procedure  (see  Section  3  for  details).  Our  basic  strategy  is  to  solve 
the  problem  by  solving  its  dual.  The  dual  problem  is  formulated  in  such 
a  way  that  it  is  an  unconstrained  maximization,  with  respect  to  the 
Lagrange  multipliers,  of  a  differentiable  concave  function  with  known, 
explicit  derivatives.  The  function  to  be  maximized  in  the  dual  itself 


is  given  as  a  certain  constrained  minimization  problem  (called  the 
augmented  primal  problem).  The  augmented  primal  problem  is  the  same 
as  the  original  (or  primal)  problem  except  that  some  of  the  constraints 
are  incorporated  to  the  objective  function  through  the  Lagrange 
multipliers.  The  dual  problem  is  solved  by  solving  the  augmented 
primal  problem  repeatedly.  This  augmented  primal  problem  is  decomposed, 
in  two  stages,  into  eventually  a  collection  of  small  size  problems  each 
of  which  corresponds  to  a  single  finite  element  under  a  single  loading 
condition. 

In  the  first  stage  of  the  decomposition,  the  augmented  primal  prob¬ 
lem  is  decomposed  into  the  "sun"  ofjjiinimization  problems  each  of  which 
is  associated  with  a  single  finite  element  under  all  the  loading 
conditions.  In  the  second  stage,  each  of  these  problems  is  partitioned 
in  such  a  way  that  the  minimization  becomes  "nested";  i.e.  the  problem 
is  transformed  into  an  "outer  minimization"  inside  of  which  many  "inner 
minimizations"  are  performed.  The  outer  minimization  is  (virtually)  an 
unconstrained  problem  with  respect  to  the  design  variables  associated 
with  a  single  finite  element,  while  each  of  the  inner  minimizations  is 
a  constrained  minimization  problem  corresponding  to  a  single  finite 
element  under  a  single  loading  condition  given  f i xed  values  of  the 
design  variables  associated  with  the  element.  We  might  note  that  the 
dual  problem  is  solved  by  means  of  a  certain  iterative  scheme  developed 
recently  ( [1 ]-[2] ) .  Figure  1  gives  a  conceptual  outline  of  our 
decomposition  procedure. 

The  present  work  is  motivated  and  inspired  by  the  following  two 
papers;  Thierauf  [3'  and  Woo  and  Schmit  [4],  Thierauf  proposes  an 


f 
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iterative  scheme  for  perfectly  plastic  optimal  design  problem  having  a 
linear  objective  function  and  quadratic  yield  conditions.  His  procedure 
has  a  flavor  of  decomposition  in  the  sense  that  in  each  iteration  an 
elastic  problem  is  solved  with  respect  to  the  finite  elements  but 
a  single  loading  condition.  For  a  general  case  with  nonlinear  objective 
and  yield  conditions,  Thierauf  suggests  repeated  linear/quadratic 
approximations  in  an  iterative  way. 

The  decomposition  procedure  proposed  by  Woo  and  Schmit  is  based  on 
a  principle  entirely  different  from  ours,  but  similar  to  ours  in  the 
sense  that  the  problem  is  decomposed  down  to  a  collection  of  many  small 
size  problems  corresponding  to  a  single  finite  element  under  a  single 
loading  condition.  The  procedure  proposed  in  this  present  paper  seems 
to  have  several  advantages  over  that  by  Woo  and  Schmit  as  explained 
below. 

The  procedure  by  Woo  and  Schmit  is  an  application  of  Dantzig-Wolf 
decomposition  principle  (Dantzig  [5])  for  (generalized)  linear  program¬ 
ming.  The  linearity  assumption  of  the  cost  function,  therefore,  is 
crucial  in  their  method  while  our  procedure  only  requires  convexity  and 
(some)  separability.  It  is  well-known  ([6])  that  the  performance  of 
algorithms  based  on  Dantzig-Wolf  decomposition  principle  is  highly 
erratic.  Further,  the  convergence  of  Woo-Schmit  procedure  is  not 
guaranteed  unless  the  yield  function  at  each  element  is  piecewise 
linear.  Our  procedure  is  guaranteed  to  converge  under  a  general 
convex  yield  function. 
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Of  course,  this  doesn't  detract  the  importance  of  the  pioneering 
work  by  Woo  and  Schmit.  In  this  early  stage  of  the  obviously  important 
research  on  the  decomposition  of  large  scale  structional  analysis,  all 
valid  approaches  must  be  tried  and  tested;  and  our  procedure  provides 
one  alternative.  Another  promising  approach  based  on  Winkler's 
factorization/decomposition  scheme  will  be  proposed  elsewhere  by  the 
first  author  of  the  present  paper  in  the  near  future. 

The  organization  of  the  rest  of  this  paper  is  as  follows.  In  the 
next  section  we  shall  specify,  as  a  mathematical  programming  problem, 
the  structural  design  problem  we  shall  consider  in  this  paper.  A 
description  of  the  proposed  decomposition  procedure  is  presented  in 
Section  3.  In  the  fourth  and  final  section  we  shall  discuss  some 
significant  simplifications  of  the  proposed  procedure  when  it  is  applied 
to  a  class  of  certain  specific  problems;  a  numerical  example  and  the 
result  of  computation  will  also  be  given. 

2 .  The  Problem 

The  problem  we  shall  be  concerned  with  can  be  stated  as  that  of 
finding  a  minimal  cost  design  for  a  perfectly  plastic  structure  sub¬ 
jected  to  several  alternative  load  conditions  under  the  overall  static 
equilibrium  condition  and  the  requirement  that  no  yieldings  occur. 

’  The  corresponding  mathematical  program  looks  like: 

I 

t 
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(1.1) 


minimize  w^y^) 


subject  to 

GjlF^.J,)  <  0 


1*1 ». . .  ,r 


y,-  >  o 


n4f 


i  i 


(1.2) 

(1.3) 

(1.4) 


Here,  is  the  vector  of  design  variables  associated  with  the  i-th 
finite  element,  PJ  is  the  vector  representing  the  j-th  loading 
condition,  is  the  vector  of  (generalized)  stresses  at  the  i-th 
element  under  the  j-th  loading  condition,  (•)  is  the  cost  of 
design  ,  (1.2)  represents  the  yielding  conditions  and  (1.4)  denotes 

the  overall  static  equilibrium. 

The  framework  (1),  with  the  assumptions  to  be  specified  later  in 
this  section,  provides  an  idealized  and  simplified,  yet  fairly  broad 
perfectly  plastic  optimization  model.  We  note  that  (1)  is  in  almost 
exactly  the  same  form  as  the  formulations  given  in  Thierauf  [3]  and 
Pape  and  Thierauf  [7],  except  for  some  generalities  allowed  in  (1). 

In  the  following,  we  shall  explain  exactly  what  these  symbols  in  (1) 
mean  and  spell  out  the  set  of  assumptions  we  shall  make.  Our  presenta¬ 
tion  emphasizes  idealized  mathematical  properties  and  requirements;  the 
(important)  question  of  how  valid  and  useful  the  formulation  (1)  is  as 
a  structural  design  model  is  left  to  expositions  by  specialists  of 
structures  (such  as  [3],  [7],  Pape  [8]  and  Rozvany  [9]).  In  particular, 
we  shall  be  tacitly  making  some  simpliefied  assumptions  such  as 
quasi  static  loading  and  "fixed  layout". 
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Consider  a  perfectly  plastic  structure  which  is  discretized  into  r 
finite  elements.  The  shape  and  nature  of  the  elements  must  be  chosen 
depending  on  the  particular  problems  and  are  irrelevant  to  our  formulation. 
Let  F.j  and  y^  denote,  respectively,  the  s..  dimensional  vector  of 
(generalized)  stresses  and  the  t^  dimensional  nonnegative  vector  of  design 
variables  associated  with  the  i-th  element,  i=l,...,r.  We  assume  that  for 
element  i  there  are  p^  yield  conditions  each  of  which  is  determined  by 
a  differentiable,  strictly  convex  yield  function,  i=l,...,r.  Thus,  the 
combined  yield  condition  for  the  i-th  element  can  be  written  as 

Gi(Fi,y.)  <  0,  1  =  1,..., r,  (2) 

where  is  a  vector  function  with  p..  components  such  that  each 
component  function  is  differentiable  and  strictly  convex  (jointly)  in 
Fj  and  y.j .  An  example  of  G..  with  p^  =  1,  t..  =  1  is 

Gi(Fi>yi>  ■  W  -  v 

where  <|>.  is  a  real -valued  "plastic  potential"  function.  Figure  2 
depicts  the  yield  condition  (2)  where  p..  =  2,  t.  =  1  and  s.  =  2. 


Fig.  2  Piecewise  differentiable  yield  function 
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The  stress  vectors  Fj,...,Fr  for  all  the  elements  are  required 
to  satisfy  the  equilibrium  condition: 


I 


r 

i=l 


NiFi  *  P’ 


(3) 


where  is  an  m  by  s^  matrix,  i=l,...,r,  and  P  is  an 
m- vector  of  external  nodal  loads.  We  assume  that  the  matrix 


N  =  (NrN2,...,Nr) 


has  full  row  rank. 

We  let  w(y1,...,yr)  denote  the  cost  (weight,  price  of  material, 
etc.)  of  design  and  assume  that  w(*)  is  separable  in  y. ,  i.e., 

w(yr...,yn)  =  w^y.),  (4) 

where  for  each  i  e  (1 . r},  w .(•)  is  a  real-valued,  differentiable 

convex  function.  Some  degree  of  non-separability  can  be  dealt  with  in 

our  model  but  we  don't  elaborate  here.  In  most  practical  cases,  however, 

the  separability  seems  a  reasonable  assumption  anyway. 

Following  Theirauf  [3]  we  formulate  the  model  so  that  it  allows 

alternative  loading  conditions.  As  Thierauf  points  out,  one  can  design 

a  structure  which  withstands  any  convex  combination  of  given  k  loading 

conditions  by  requiring  k  sets  of  equilibrium  equations.  To  accomo- 

1  2  k 

date  this,  we  let  m-vectors  P  ,P  ,...,P  represent  the  k  loading 
conditions,  respectively,  and  require  (in  place  of  (3)): 


N,F 


i’  i 


1-1,. 


»r 


(5) 
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for  each  j=l,...,k,  where  for  each  i  <■  (1 , . . .  ,rl ,  F'j  is  the  s^-vector 
of  stresses  associated  with  element  i  under  the  j-th  load  condition. 
Accordingly,  we  replace  the  yield  condition  (1)  with 

Gi(F'j,yi)  <  0,  j=l,...,k, 

for  each  i  c  {1 . r}.  The  overall  optimization  problem  is  thus  given 

by  (1).  For  notational  simplicity  we  shall  use  the  symbols  y  and  F, 
respectively,  to  denote  the  super  vectors. 


y  =  (yr-.*»yr)  and 
F  =  (F^  F^  F*"  F^ 


Finally,  let  us  note  the  size  of  (1).  The  problem  has 
IjTl  (t-j+  k*sj)  variables,  nonlinear  and  k*n  linear 

constraints  (plus  jT,  t^  nonnegativity  conditions).  Recall  that: 


r  =  #  of  finite  elements 

m  =  #  of  overall  equilibrium  equations 

k  =  r  of  loading  conditions 

t.  =  #  of  design  variables  for  element  i 

s.j  =  #  of  stress  componenets  for  element  i ;  and 

p.  =  *  of  yield  conditions  for  element  i. 


For  every  i,  each  of  t^ ,  s^  and  p.  is  relatively  small.  Also,  k 
is  expected  to  be  relatively  small.  The  number  of  elements,  r,  is  by 
far  the  most  crucial  quantity  in  a  large  scale  problem;  it  could  be 
hundr  Js  or  even  thousands. 
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overcome  this  difficulty,  we  shall  employ  the  proximal  point  algorithm 
(Rockafellar  [1]);  i.e.,  we  augment  the  objective  function  by  adding 
certain  quadratic  terms  to  induce  "nice"  properties  in  the  augmented 
primal  problem,  such  as  the  uniqueness  of  the  solution  and  the  dif¬ 
ferentiability  of  h  (see  below),  and  solve  the  augmented  problem 
iteratively,  until  the  true  solution  is  found.  Obviously,  solving 
the  problem  iteratively  is  costly,  but  the  benefit  of  having  the  nice 
properties  far  outweighs  the  added  complication. 

More  specifically,  for  a  given  vector  (y,F),  we  consider  the 
problem  (1)  where  the  objective  function  (1.1)  is  replaced  with 


where  X  is  some  given  positive  scalar  and  IMI  is  the  Euclidean 
2-norm.  Thus,  instead  of  the  dual  problem  (6)-(7),  we  shall  consider 

maximize  h(u)  (8) 

u  =(u1,...,uk) 

with 


>  o 


VFi’*i>  £° 

i=l,...,r,j=l,...,k 


t 

> 
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where  we  set 


ff(yf)  =  +  2X  lly^y^l2,  and  (10.1) 


j-1 . . ,k. 

It  must  be  noted  that  f..  and  g^  are  functions  of  y^  and  F'j 
as  well  as  y^  and  F1?,  respectively,  and  that  in  fact  y  and  F 
are  changed  from  one  iteration  to  another;  but  we  suppress  them  for 
notational  simplicity. 

Let  u*  be  the  dual  optimum  and  let  (y*,F*)  be  the  correspond¬ 
ing  minimizer  in  (9).  If  (y*,F*)  coincides  with  (y,F),  then  it 
can  be  shown  (compare  (6)-(7)  and  (8)-(9))  that  (y*,F*)  is  the 
optimal  solution  to  the  original  problem  (1).+  Otherwise,  we  replace 
(y,F)  with  (y*,F*)  and  repeat  the  process.  This  process  is 
guaranteed  to  converge,  with  a  reasonably  good  convergence  rate^ 

(see  [2])  and  computational  results  (see  Ha  [10])  have  been  encourag¬ 
ing.  This  iterative  process  is  the  basic  iteration  scheme  providing 
the  "highest  echelon"  loop  in  our  decomposition  procedure.  The 
augmented  primal  problem  in  (9)  will  be  decomposed,  in  two  stages, 
into  eventually  small  size  problems  each  of  which  is  associated  with 
a  single  structural  element  under  a  single  loading  condition;  but  we 
shall  first  examine  the  properties  of  the  augmented  primal  further. 


4.  + 

Of  course,  in  practice,  we  accept  a  solution  (y*,F  )  if  it  is 
"sufficiently"  close  to  (y,F). 

4*4. 

If  w..(y^)  is  linear  for  each  i,  then  the  process  terminates 
after  a  finite  number  of  iterations. 
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For  each  u,  let  (y(u),  F(u))  denote  the  unique  solution  of  the 
augmented  primal  problem.  It  can  be  shown  (see  [11])  that  n(u)  is 
strictly  concave  and  differentiable  in  u.  Furthermore,  the  derivatives 
of  h(u)  is  known  explicitly;  i.e.,  the  gradient  vector  of  h(u)  with 
respect  to  u^  is  given  by 

vuj  Mu)  =  I ^1  N.F^u)  -  PJ, 

where  F'j(u)  is  the  appropriate  component  of  F(u)  defined  above.  We 
note  that  all  these  crucial  and  nice  properties  are  the  results  of  adding 
the  quadratic  terms.  The  maximization  problem  (8)  does  have  a  large  size 
(the  dimension  of  u  is  k*m);  however,  the  facts  that  h(u)  is 
differentiable,  strictly  concave  with  known,  explicit  derivatives  and 
that  the  problem  is  unconstrained  ensure  that  the  optimization  can  be 
performed  effectively  by  using  some  of  the  most  reliable  and  fast  uncon¬ 
strained  nonlinear  programming  routines  (such  as  a  quasi -Newton  method). 

The  first  stage  of  the  decomposition  process  is  done  by  observing 
that  the  optimization  in  (8)  can  be  accomplished  by  solving  the  following 
subproblem  separately  for  each  i: 

minimize  f.^)  +  J  j  (^  )  +  (uJ )  N^)  (11.1) 

j  ^ 

subject  to  (F:j,y.j )  <  0,  j=l . k  (11.2) 


y,  >  o- 


(11.3) 
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This  is  easy  to  see  if  we  rearrange  some  terms  in  (9)  and  rewrite  it  as 


h(u)  = 


r  k  •  -i-  j 

V  V  rJ 


Vi  >  o 

Gi(Fj,yi)  <  0 
i=l , . . . ,r,j=l . .  ,k 


{  KMy.)*  l  (gn(H))+  l  1  (uj)'n.fJ 

i=l  1  1  j=l  1J  1  i=l  j«l  1  1 


-  I(UJ)V) 

j=l 


■  -  2  (uj>tpj  *  2 

j=l  i= 


p  k  •  •  • 

.1  min  {f  1  (yn- )+  M9ij-(F^)+(uj)TNiF^)}. 

yi  >  o 
61(F^^i)  <  0 
j=l .  ,k 


Note  that  the  problem  (11)  corresponds  to  a  single  finite  element  and 
all  k  loading  conditions. 

The  second  stage  of  the  decomposition  is  obtained  by  re-expressing 
the  problem  (11)  as  the  following  nested  minimization  (it  is  not 
difficult  to  verify  the  validity  of  this  "partitioning"): 


minimize  (f,(y.)+  l  v..(yj} 
y,.  >  0  11  j=l  13  1 


where 


vij(yi)  =  min  { g. j )+( )TN. F^ > .  ( 
G.(F3,yi)<0 

Namely,  (11)  is  solved  by  performing  the  outer  minimization  (12)  with 
respect  to  y^ ,  where  the  value  of  each  v^.(y^)  is  given  by  solving 
the  following  inner  ninimization: 


k  # 


minimize 


Sijd’)  +  (uj)th1fJ 


(14.1) 


subject  to  G.(F^,y. )  <  0.  (14.2) 

It  can  be  shown  (see  Rockafellar  [11])  that  the  function 

f i  (yi )  +  ^j=l  vij^yi^ 

is  convex  in  y. .  Further,  under  certain  mild  conditions  (see  Ha  [10]) 
this  function  is  differentiable  in  y^  and  an  explicit  form  its 
derivatives  can  be  obtained.  These,  as  well  as  the  fact  that  the 
dimension  of  y^  is  small,  lead  to  an  efficient  minimization  (using 
a  quasi-Newton  algorithm,  for  example). 

Now,  the  inner  minimization,  (14),  is  a  (relatively)  small  size 
constrained  optimization  problem  with  s^  variables  and  p^  con¬ 
straints,  which  corresponds  to  a  single  finite  element  under  a  single 
loading  condition  for  given  values  of  the  design  variables.  By 
hypothesis,  G-(F^.y^)  is  differentiable  and  convex  in  F^  (for 
fixed  y.j);  the  objective  function  is  quadratic  (c.f.  (10)).  These 
properties  and  the  smallness  of  the  problem  ensures  an  efficient 
solution  of  the  inner  minimization. 

To  summarize  the  proposed  procedure,  we  look  at  it  as  a  "modular 
computer  package"  consisting  of  several  subroutines  as  follows. 

Preparation: 

Choose  an  initial  point  (y,F).  Add  the  quadratic  terms  to  the 
objective  function  (as  in  (10)). 

I 
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Main  Routine  (Basic  Iteration  Scheme) 

Call  subroutine  MAXH  to  solve  the  problem  (8)-(9).  Let  (y*,F*)  be 
the  optimal  solution  to  the  augmented  primal  problem  (9)  corresponding 
to  the  optimal  solution  u*  of  (8).  If  ll(y*,F*)-  (y,F)||  <  e  for  a 
prescribed  tolerance  e  >  0,  then  stop  ((y*»F*)  is  optimal); 
otherwise,  set  y  =  :  y*,  F  *  :  F*  and  repeat. 

Subroutine  MAXH  (Input:  (y,F)  Output:  (y*,F*)) 

Use  a  quasi-Newton  algorithm  to  solve  (8)-(9);  call  subroutine  EVAHDH 
to  evaluate  h(u)  and  its  derivatives  needed  for  the  algorithm. 

Subroutine  EVAHDH  (Input:  (y,F),  u  Output:  h(u),  Vh(u)) 

Call  subroutine  OUTMIN  to  solve  (12)-(13)  for  each  i=l,...,r;  let 
12  w 

(y.j(u);  F..(u),  F^(u),. . . ,F^(u))  be  the  solution,  i=l,...,r.  The 
values  of  h(u)  and  its  derivatives  are  given,  respectively,  by: 


V  i  h(u)  «  l  N.FJ(u)  -  Pj,  j=l ,. . . ,k. 
uJ  i=l  1  1 

Subroutine  OUTMIN  (Input:  (y,F),  u,  i  Output:  (y(u),  F(u ) ) ) 

Solve  the  problem  (12)  for  y^ ;  call  subroutine  INMIN  to  evaluate 


v,,(yJ,  j=l . k. 


Subroutine  INMIN  (Input:  F,  u,  i,  y^,  j  Output:  v - j (y^ )) 
Solve  the  problem  (14). 

Figure  3  shows  the  relationships  among  these  subroutines. 


Basic  Iteration  Scheme 


The  decomposition  procedure  was  presented  in  the  previous  section 


for  a  fairly  broad  class  of  problems.  Most  often,  the  particular 
problem  under  consideration  has  many  special  features  which  lead  to 
some  simplifications  of  the  procedure.  In  this  section  we  shall 
examine  one  such  case. 

The  problem  we  shall  consider  here  is  given  by: 


minimize 
subject  to 


?<Ft>  Vi  S  *1 


y i  >  o 

I  =  Pj 
i=1  1  1 


i=l 

j=l 


,r 

»k 


(15.1) 

(15.2) 

(15.3) 

(15.4) 


where  D..  is  an  s..  by  s.  symmetric  positive  definite  matrix  for 
each  ie{l,...,r}.  This  is  the  special  case  of  the  general  problem 
(1),  where  ti  =  Pi  =  1  for  each  i  and  is  the  real-valued 
function  with 

Gi(Fi‘  =  }(Fi)TDiFi  -  y,’  1*1  »•  •  •  »r. 

Namely,  each  finite  element  has  a  single  design  variable  and  a  single, 
quadratic  yield  condition.  In  this  case,  rather  significant  simplifi¬ 
cations  can  be  made  in  the  proposed  decomposition  procedure  particu¬ 
larly  with  respect  to  solving  (12)-(13). 
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Consider  the  nested  minimization  problem  (12)-(13)  (here,  u  is 
given).  Since  y^  is  a  scalar,  the  outer  minimization  is  a  one¬ 
dimensional  optimization,  or  a  "line  search".  It  can  be  shown  (not 
difficult)  that  the  function  to  be  minimized  is  differentiable  and 
convex  where  the  derivative  is  equal  to  the  value  of  optimal  Lagrange 
multiplier.  Thus,  the  outer  minimization  can  be  performed  extremely 
efficiently  by  using  a  one-dimensional  quasi-Newton  algorithm.  By  using 
other  properties  of  v..(»)»  further  simplifications  are  possible  (but 

*  J 

not  elaborated  here).  It  turns  out,  also,  that  the  inner  minimization 
to  evaluate  v.^y.)  reduces  to  solving  one  (nonlinear)  equation  with 
one  unknown  and  so  it  can  be  performed  in  an  extremely  simple  manner. 
Details  are  described  at  the  end  of  this  section. 

We  wrote  a  Fortran  computer  code  to  solve  an  optimal  design 
problem  of  the  type  (15)  using  the  proposed  decomposition  procedure  and 
the  simplifications  in  solving  (12)- (13)  outlined  above.  We  applied  the 
algorithm  to  solve  a  numerical  example  treated  in  Thierauf  [3] 

(Example  1  on  p.  146)  and  it  computed  an  optimal  solution  after  4  "outer" 
iterations(Basic  Iteration  Scheme)  and  38  quasi-Newton  steps  (75  function 
evaluations).  The  solution  is  feasible  up  to  the  error  of  magnitude  10"4 
and  its  objective  value  is  slightly  lower  than  that  given  in  Thierauf's 
paper.  We  have  solved  a  few  other,  larger  size  problems  using  the  algo¬ 
rithm  also. 

In  the  remainder  of  this  section  we  shall  explain  how  to  solve 
the  problem  (14)  efficiently  when  the  problem  is  given  by  (15).  For 
notational  simplicity  we  shall  write  the  problem  as  follows: 

minimize  ^oJx  -  c^x  (16.1) 

subject  to  ^-xTDx  <  a,  (16.2) 


jc 


where  x  is  an  s-vector  variable,  c  is  a  constant  s-vector,  D  is 
a  constant  symmetric  and  positive  define  matrix  and  a  is  a  given 
positive  scalar.  Clearly  the  inner  minimization  problem  (14)  for  (15) 
is  obtained  by  setting 

x  =  F? ,  D  =  Di ,  a  =  yi , 

c  =  -  <u1>V 

and  by  adding  the  constant  ^  (F^)TF'j  to  the  objective  of  (16). 

The  problem  is  obviously  feasible  and  is  a  convex  program  with  a 
strictly  convex  objective  function;  thus  it  has  a  unique  minimum, 


which  is  characterized  by  the  Kuhn-Tucker  conditions: 


x  -  c  +  £Dx  =  0 


C  >  0,  x  Dx  <  a 


£  •  (a-x  Dx)  =  0, 


(17.1) 

(17.2) 


(17.3) 


where  £  is  the  scalar  Lagrange  multiplier.  Clearly,  the  optimal 
Lagrange  multiplier  is  zero  if  and  only  if  cTDc  <  a;  if  £  =  0, 
then  x  =  c  is  the  optimal  solution  to  (16).  Thus,  the  first  step 
of  the  procedure  to  solve  (16)  should  be  to  check  if  c^Dc  <  a 
holds  (if  so  x  =  c  is  the  optimum). 

Suppose,  now,  that  x  *  c  is  not  feasible  in  (16)  and  so  the 
optimal  Lagrange  multiplier  is  strictly  positive.  To  solve  (16) 
under  this  circumstance,  we  take  the  dual  approach  as  we  did  for  the 
original  problem  (1);  the  dual  of  (16)  is  given  by 
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maximize  a(C),  (18) 

C  >  0 


where 


a(0 

For  each  fixed  C> 
it  is  not  difficult 

x(0 

where  I  is  the  s 

a  (O 


=  min  {IoJx-c^x+icx^Dx}  -  Ca. 
x 

let  x(C)  denote  the  minimizer  in  (19). 
to  see  that 

•  (I+CD)"^, 

by  s  identity.  Hence  we  have 
=  ]*(C)Tx(C)  -  cTx(C)  +  Cx(C)TDx(C)  -  Cot 
=  ^-x(c)T[x(c)  -  c  +  CDx(C)]  -  ^cTx(C)  -  Cot 


(19) 


Then 


(20) 


=  -  |-cTx  ( C )  -  Ca 
=  -  ]-cT(I+C0)_1c  -  Ca. 


As  before,  a(C)  is  differentiable  concave  in  C  and  so  the 
dual  problem  (18)  can  be  solved  by  obtaining  a  nonnegative  (actually 
positive)  solution  C  to  the  single  equation: 

^  []<T(I+CD)'1c+Ca]  =  0, 

or 

^CT(I+CD)_1C  =  -2a.  (21) 


xiAiliiJO KaZTT' 
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For  instance,  if  s  =  1,  then  the  equation  (21)  becomes 

2 

d  zc  _ 

TThdT  '  "2“ 

where  both  c  and  D  are  scalars  and  thus  the  positive  root  of 
equation  is  given  by+ 


€  - 


In  a  general  case  an  alternative  form  of  the  equation  (21)  may 
be  obtained  as  follows.  Let  {t^  : i = 1 . .  ,s}  and  (h1 :i=l , . . . ,s}  be 
the  sets  of  eigenvalues  and  eigenvectors,  respectively,  of  the  matrix 
D.  By  the  assumptions  on  D,  r^.  is  a  positive  scalar  and  h1  is 
real,  for  each  i.  It  is  not  difficult  to  show  that  h1  is  also  an 
eigenvector  of  the  matric  (I+£D)”^  with  the  associated  eigenvalue 
l/(l+£;n.)  for  each  i.  Let  y  =  (y1,....y$)T  be  the  unique  solution 
of  the  system 

c  *  i  V«ih1 ;  or 

y  =  [h1 , . . . ,hs]_1c.  (22) 

Then  it  can  be  shown  (not  difficult)  that 


xU) 


=  (I+SD)-1c 


S  M-.  j 


- -  ^ 

'  This  E,  is  positive  if  and  only  if  c  0/2  >  a;  if  c  D/2  <  a,  or 

cDc/2  <  a,  then  x  =  c  is  feasible  and  thus  the  problem  would  be 

solved  by  the  "first  step",  i.e.  the  check  of  feasibility  of  x  =  c. 
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Therefore,  the  equation  (21)  becomes 


I  (cV) 
=  1 


(1+Cni )2 


=  2a. 


This  process  requires  the  determination  of  all  eigenvalues  and  eigen¬ 
vectors  of  D  and  the  matrix  inversion  necessary  in  (22).  Most  comput¬ 
ing  centers  are  equipped  with  computer  packages  to  do  these  operations. 
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